proj cub fit

import%20marimo%0A%0A__generated_with%20%3D%20%220.10.8%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20plotly.graph_objects%20as%20go%0A%20%20%20%20import%20scipy%20as%20scp%0A%20%20%20%20import%20itertools%20as%20itt%0A%20%20%20%20from%20typing%20import%20NamedTuple%0A%20%20%20%20return%20NamedTuple%2C%20go%2C%20itt%2C%20mo%2C%20np%2C%20scp%0A%0A%0A%40app.cell%0Adef%20_(NamedTuple%2C%20np)%3A%0A%20%20%20%20class%20AffineMap(NamedTuple)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Represents%20an%20affine%20transformation%20x%20%5Cmapsto%20Ax%20%2B%20b%20on%20R%5E(n)%0A%20%20%20%20%20%20%20%20via%20its%20associated%20linear%20transformation%20A%20and%20offset%20b.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%0A%20%20%20%20%20%20%20%20linear_transf%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20offset%3A%20np.ndarray%0A%0A%20%20%20%20%20%20%20%20def%20apply_to(self%2C%20points)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20self.linear_transf%20%40%20points%20%2B%20self.offset.reshape(3%2C%201)%0A%0A%20%20%20%20%20%20%20%20def%20__matmul__(self%2C%20other)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20type(other)%20is%20AffineMap%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20self.after(other)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20self.apply_to(other)%0A%0A%20%20%20%20%20%20%20%20def%20after(self%2C%20other%3A%20%22AffineMap%22)%20-%3E%20%22AffineMap%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22Computes%20the%20composition%20self%20%5Ccirc%20other.%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20AffineMap(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.linear_transf%20%40%20other.linear_transf%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.linear_transf%20%40%20other.offset%20%2B%20self.offset%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20def%20inv(self)%20-%3E%20%22AffineMap%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20inv%20%3D%20np.linalg.inv(self.linear_transf)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20AffineMap(inv%2C%20-inv%20%40%20self.offset)%0A%20%20%20%20return%20(AffineMap%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Some%20helpers%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20Let%20%24v_1%2C%20...%2C%20v_k%24%20be%20%24k%24%20linearly%20independent%20vectors%20in%20%24%5CR%5En%24.%0A%20%20%20%20%20%20%20%20We%20want%20to%20determine%20a%20matrix%20%24A%24%20such%20that%20%24A%20v_i%20%3D%200%24%20for%20all%20%24i%24%20and%20%24A%20x%20%3D%20x%24%20for%20all%20%24x%20%5Cin%20W%5E%5Cbot%24%20where%20%24W%20%3A%3D%20%5Coperatorname%7Bspan%7D(v_1%2C%20...%2C%20v_k)%24.%0A%0A%20%20%20%20%20%20%20%20Let%20%24V%20%3D%20(v_1%2C%20...%2C%20v_k)%20%5Cin%20%5CR%5E%7Bn%2Ck%7D%24.%0A%20%20%20%20%20%20%20%20It's%20a%20standard%20result%20that%20%24P%20%3A%3D%20V%20V%5E%5Cdagger%24%20is%20the%20orthogonal%20projector%20onto%20%24W%24%20---%20here%20%24V%5E%5Cdagger%24%20denotes%20the%20moore-penrose%20pseudoinverse%20of%20%24V%24.%20Since%20the%20%24%5C%7Bv_i%5C%7D%24%20are%20linearly%20independent%20%24V%24%20has%20full%20column-rank%20and%20hence%20we%20have%20%24V%5E%5Cdagger%20%3D%20(V%5ET%20V)%5E%7B-1%7D%20V%5ET%24%20and%20hence%20%24P%20%3D%20V%20(V%5ET%20V)%5E%7B-1%7D%20V%5ET.%24%0A%0A%20%20%20%20%20%20%20%20Now%20let%20%24V%20%3D%20U%20%5CSigma%20Q%5ET%24%20be%20a%20SVD%20of%20%24V%24.%0A%0A%20%20%20%20%20%20%20%20Then%20%24V%20(V%5ET%20V)%5E%7B-1%7D%20V%5ET%20%3D%20V%20(Q%20%5CSigma%5ET%20U%5ET%20U%20%5CSigma%20Q%5ET)%5E%7B-1%7D%20V%5ET%20%3D%20V%20(Q(%5CSigma%5ET%20%5CSigma)%20Q%5ET)%5E%7B-1%7D%20V%5ET%20%3D%20V%20Q%20(%5CSigma%5ET%20%5CSigma)%5E%7B-1%7D%20Q%5ET%20V%5ET%20%3D%20U%20%5CSigma%20(%5CSigma%5ET%20%5CSigma)%5E%7B-1%7D%20%5CSigma%5ET%20U%5ET%24.%20Since%20%24%5CSigma%5ET%20%5CSigma%24%20is%20a%20diagonal%20matrix%20it%20follows%20that%20%24P%20%3D%20UU%5ET%24%20and%20hence%20%24A%20%3A%3D%20I%20-%20UU%5ET%24%20works.%0A%0A%20%20%20%20%20%20%20%20Similarly%20we%20can%20show%20that%20for%20a%20thin%20%2F%20reduced%20QR%20decomposition%20%24V%20%3D%20QR%24%20we%20have%20%24A%20%3D%20I%20-%20QQ%5ET%24.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20proj_onto_orthogonal_complement_single(v)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Similar%20to%20projon%E2%86%92orthogonalcomp%E2%89%A4mentproj_onto_orthogonal_complement%20special-cased%20to%20a%20single%20vector.%0A%0A%20%20%20%20%20%20%20%20Note%3A%20One%20could%20also%20implement%20a%20special%20version%20for%20two%20vectors%20in%20R%5E3%20or%20more%0A%20%20%20%20%20%20%20%20generally%20n-1%20vectors%20in%20R%5En.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20v_unit%20%3D%20v%20%2F%20np.linalg.norm(v)%0A%20%20%20%20%20%20%20%20projection_matrix%20%3D%20np.outer(v_unit%2C%20v_unit)%0A%20%20%20%20%20%20%20%20return%20np.eye(v.shape%5B0%5D)%20-%20projection_matrix%0A%0A%0A%20%20%20%20def%20mat_from_eigenpairs(*eigenpairs%3A%20tuple%5Bfloat%2C%20np.ndarray%5D)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Computes%20a%20matrix%20A%20such%20that%20Av%20%3D%20tv%20for%20a%20given%20collection%20of%20pairs%20(t%2Cv).%0A%0A%20%20%20%20%20%20%20%20This%20solves%20systems%20Q%5E(T)%20a_i%5E(T)%20%3D%20Lambda%20Q%5E(T)%20where%20Q%20and%20Lambda%20are%20matrices%20as%20in%20the%20eigendecomposition%20and%20a_i%20is%20the%20i-th%20row%20of%20A.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20eigenvals%20%3D%20np.array(%5Bval%20for%20(val%2C%20vec)%20in%20eigenpairs%5D)%0A%20%20%20%20%20%20%20%20eigenvecs%20%3D%20np.row_stack(%5Bvec%20for%20(val%2C%20vec)%20in%20eigenpairs%5D)%0A%20%20%20%20%20%20%20%20if%20eigenvecs.shape%5B0%5D%20!%3D%20eigenvecs.shape%5B1%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Need%20n%20eigenpairs%20to%20determine%20matrix%20acting%20on%20R%5En.%22%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20scaled%20%3D%20eigenvals.reshape(-1%2C%201)%20*%20eigenvecs%0A%20%20%20%20%20%20%20%20mat%20%3D%20np.row_stack(%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.linalg.solve(eigenvecs%2C%20scaled%5B%3A%2C%20i%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(scaled.shape%5B1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20mat%0A%0A%0A%20%20%20%20def%20proj_onto_orthogonal_complement(*vectors%3A%20np.ndarray)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Returns%20a%20linear%20transformation%20that%20projects%20onto%20the%20orthogonal%20complement%20of%20the%0A%20%20%20%20%20%20%20%20span%20of%20a%20given%20collection%20of%20vectors%20%5C%7Bv_i%5C%7D.%20The%20vectors%20are%20assumed%20to%20be%20linearly%0A%20%20%20%20%20%20%20%20independent.%0A%0A%20%20%20%20%20%20%20%20So%20it%20maps%20all%20the%20vectors%20v_i%20to%200%20and%20(span%7Bv_i%7D_i)%5E%5Cbot%20is%20invariant%20under%20the%20map.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20%23%20we'd%20really%20want%20to%20return%20the%20identity%20if%20len(vectors)%20%3D%3D%200%20but%20in%20this%20case%20we%20don't%20know%20the%20dimension%0A%20%20%20%20%20%20%20%20if%20len(vectors)%20%3D%3D%201%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20this%20isn't%20really%20necessary%20and%20the%20below%20code%20works%20for%20a%20single%20vector%20but%20eh%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20proj_onto_orthogonal_complement_single(vectors%5B0%5D)%0A%20%20%20%20%20%20%20%20%23%20if%20len(vectors)%20%3D%3D%202%3A%20%23%20for%203D%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20return%20mat_from_eigenpairs(%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20%20%20%20%20(0%2C%20vectors%5B0%5D)%2C%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20%20%20%20%20(0%2C%20vectors%5B1%5D)%2C%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20%20%20%20%20(1%2C%20np.cross(vectors%5B0%5D%2C%20vectors%5B1%5D))%2C%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20vec_mat%20%3D%20np.column_stack(vectors)%0A%20%20%20%20%20%20%20%20%20%20%20%20qr_decomp%20%3D%20np.linalg.qr(vec_mat%2C%20mode%3D%22reduced%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20projection_matrix%20%3D%20qr_decomp.Q%20%40%20qr_decomp.Q.T%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20np.eye(vec_mat.shape%5B0%5D)%20-%20projection_matrix%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20mat_from_eigenpairs%2C%0A%20%20%20%20%20%20%20%20proj_onto_orthogonal_complement%2C%0A%20%20%20%20%20%20%20%20proj_onto_orthogonal_complement_single%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Central%20algorithm%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(AffineMap%2C%20NamedTuple%2C%20np%2C%20proj_onto_orthogonal_complement)%3A%0A%20%20%20%20class%20CubeFit(NamedTuple)%3A%0A%20%20%20%20%20%20%20%20%23%20remaps%20the%20parallelotope%20in%20R%5Ed%20such%20that%20d%20of%20its%20diagonals%20are%20orthonormal%20post-transformation%20(and%20the%20center%20is%20at%20the%20origin).%0A%20%20%20%20%20%20%20%20diagonal_orthonormalizer%3A%20AffineMap%0A%20%20%20%20%20%20%20%20%23%20maps%20the%20parallelotope%20into%20the%20unit%20cube%20%5B0%2C1%5D%5Ed.%0A%20%20%20%20%20%20%20%20cubifier%3A%20AffineMap%0A%0A%0A%20%20%20%20def%20proj_cube_fit(%0A%20%20%20%20%20%20%20%20points%3A%20np.ndarray%2C%20points_centered%3A%20bool%20%3D%20False%0A%20%20%20%20)%20-%3E%20AffineMap%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%0A%20%20%20%20%20%20%20%20Remark%3A%20for%20more%20details%20on%20the%20implementation%20please%20refer%20to%20the%20%5Bprototypes%2Fiter_proj_cube_fit_2.py%5D%20notebook.%0A%20%20%20%20%20%20%20%20The%20basic%20idea%20of%20the%20algorithm%20is%20to%20successively%20determine%20diagonals%20by%20determining%20maximal%20elements%20in%20those%20subspaces%20that%20are%20orthogonal%20to%20the%20already%20determines%20diagonals.%0A%20%20%20%20%20%20%20%20It%20utilizes%20that%20the%20image%20of%20a%20parallelotope%20under%20a%20linear%20map%20is%20again%20a%20parallelotope%2C%20that%20the%20preimages%20of%0A%20%20%20%20%20%20%20%20vertices%20in%20the%20image%20are%20also%20vertices%20in%20the%20original%20parallelotope%2C%20and%20that%20any%20maximizer%20of%20the%20norm%20over%20a%0A%20%20%20%20%20%20%20%20parallelotope%20is%20necessarily%20at%20a%20vertex%20(as%20is%20easily%20seen%20by%20noting%20that%20the%20euclidean%20norm%20is%20strongly%20convex%0A%20%20%20%20%20%20%20%20and%20hence%20f(x)%20%3D%20%7C%7CT(x)%7C%7C_2%20is%20as%20well%20for%20affine%20maps%20T%3B%20it%20also%20seems%20geometrically%20rather%20self-evident.)%0A%0A%20%20%20%20%20%20%20%20Args%3A%0A%20%20%20%20%20%20%20%20*%20points%3A%20Array%20of%20shape%20(d%2Cn)%20where%20n%20is%20the%20number%20of%20points%20and%20d%20the%20dimension%20of%20space%0A%20%20%20%20%20%20%20%20*%20points_centered%3A%20whether%20the%20inputs%20points%20may%20be%20assumed%20to%20be%20centered%20around%20the%20origin.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20dim_space%20%3D%20points.shape%5B0%5D%0A%20%20%20%20%20%20%20%20%23%20center%20points%20as%20needed%0A%20%20%20%20%20%20%20%20if%20points_centered%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20center%20%3D%20np.zeros(dim_space)%0A%20%20%20%20%20%20%20%20%20%20%20%20centered_points%20%3D%20points.copy()%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20center%20%3D%20np.mean(points%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20%20%20%20%20centered_points%20%3D%20points%20-%20center.reshape(dim_space%2C%20-1)%0A%0A%20%20%20%20%20%20%20%20%23%20determine%20d%20diagonals%20of%20the%20parallelotope%20ordered%20by%20their%20lengths%0A%20%20%20%20%20%20%20%20norms_1%20%3D%20np.linalg.norm(centered_points%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20diag_1%20%3D%20centered_points%5B%3A%2C%20np.argmax(norms_1)%5D%0A%20%20%20%20%20%20%20%20diags%20%3D%20%5Bdiag_1%5D%20%20%23%20could%20preallocate%20for%20dim_space%20vectors%20here%0A%20%20%20%20%20%20%20%20%23%20we%20want%20to%20%2F%20are%20able%20to%20determine%20n%20diagonals%20in%20total%2C%20but%20we%20already%20know%20one%0A%20%20%20%20%20%20%20%20for%20_%20in%20range(dim_space%20-%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20ortho_projector%20%3D%20proj_onto_orthogonal_complement(*diags)%0A%20%20%20%20%20%20%20%20%20%20%20%20projected_points%20%3D%20ortho_projector%20%40%20centered_points%0A%20%20%20%20%20%20%20%20%20%20%20%20norms%20%3D%20np.linalg.norm(projected_points%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20%20%20%20%20new_diag%20%3D%20centered_points%5B%3A%2C%20np.argmax(norms)%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20diags.append(new_diag)%0A%0A%20%20%20%20%20%20%20%20diags_mat%20%3D%20np.column_stack(diags)%0A%20%20%20%20%20%20%20%20%23%20This%20estimates%20the%20transformation%20by%20that%20matrix%20A%20that%20makes%20the%20diagonals%20A-orthonormal%0A%20%20%20%20%20%20%20%20%23%20and%20moreover%20maps%20the%20i-th%20diagonal%20to%20the%20i-th%20standard%20unit%20vector%0A%20%20%20%20%20%20%20%20estimated_transformation_mat%20%3D%20np.linalg.inv(diags_mat)%0A%20%20%20%20%20%20%20%20diagonal_orthonormalizer%20%3D%20AffineMap(estimated_transformation_mat%2C%20-center)%0A%20%20%20%20%20%20%20%20%23%20the%202D%20case%20would%20really%20be%20done%20at%20this%20point%20(since%20we%20can%20just%20pick%20one%20diagonal%20point%0A%20%20%20%20%20%20%20%20%23%20and%20we%20know%20that%20it's%20connected%20to%20both%20points%20of%20the%20other%20diagonal%20by%20an%20edge%3B%20in%0A%20%20%20%20%20%20%20%20%23%20constrast%20to%20the%203D%20case)%0A%20%20%20%20%20%20%20%20%23%20and%20the%20%3E3%20dimensional%20one%20needs%20a%20bit%20more%20thinking%3A%0A%20%20%20%20%20%20%20%20%23%20We%20can%20probably%20again%20pick%20the%20point%20of%20maximal%20norm%20in%20the%20transformed%20space%20as%20next%0A%20%20%20%20%20%20%20%20%23%20diagonal%20and%20then%20construct%20connected%20vertices%20similar%20to%203D%20case.%20But%20I%20haven't%20thought%20it%20out%20properly%20and%20I'm%20not%20sure%20whether%20that%20indeed%20uniquely%20determines%20the%20parallelotope%0A%0A%20%20%20%20%20%20%20%20transformed_points%20%3D%20estimated_transformation_mat%20%40%20centered_points%0A%20%20%20%20%20%20%20%20positive_first%20%3D%20transformed_points%5B0%2C%20%3A%5D%20%3E%200%0A%20%20%20%20%20%20%20%20positive_second%20%3D%20transformed_points%5B1%2C%20positive_first%5D%20%3E%200%0A%20%20%20%20%20%20%20%20positive_third%20%3D%20transformed_points%5B2%2C%20positive_first%5D%20%3E%200%0A%0A%20%20%20%20%20%20%20%20p%20%3D%20transformed_points%5B%3A%2C%20positive_first%5D%0A%20%20%20%20%20%20%20%20transformed_diag_4%20%3D%20p%5B%3A%2C%20np.argmax(np.linalg.norm(p%2C%20axis%3D0))%5D%0A%20%20%20%20%20%20%20%20%23%20for%20topological%20reasons%20the%20new%20vertex%20must%20connect%20to%20exactly%20these%20other%20vertices%0A%20%20%20%20%20%20%20%20_sgn_d4%20%3D%20np.sign(transformed_diag_4)%0A%20%20%20%20%20%20%20%20transformed_diagonals%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B1%2C%200%2C%200%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B0%2C%20_sgn_d4%5B1%5D%2C%200%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B0%2C%200%2C%20_sgn_d4%5B2%5D%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20transformed_diag_4%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20%23%20black%20magic%20or%20something%20like%20that%0A%20%20%20%20%20%20%20%20%23%20this%20just%20computes%202D%20dot%20products%20of%20various%20subsets%20of%20the%20data%20projected%20down%20onto%20a%202D%20plane.%20It%20ignores%20%22negative%20contributions%22%0A%20%20%20%20%20%20%20%20%23%20assuming%20enough%20sampling%20the%20negative%20ones%20shouldn't%20matter%20but%20I'm%20not%20quite%20sure%20how%20it%20works%20out%20if%20there's%20not%20*that*%20many%20points%20to%20begin%20with%0A%20%20%20%20%20%20%20%20%0A%20%20%20%20%20%20%20%20p%20%3D%20transformed_points%5B1%3A%2C%20positive_first%5D%0A%20%20%20%20%20%20%20%20options%20%3D%20%5B%5B1%2C%201%5D%2C%20%5B1%2C%20-1%5D%2C%20%5B-1%2C%201%5D%2C%20%5B-1%2C%20-1%5D%5D%0A%20%20%20%20%20%20%20%20p2%20%3D%20p%5B0%2C%20positive_second%5D.sum()%0A%20%20%20%20%20%20%20%20p3%20%3D%20p%5B1%2C%20positive_third%5D.sum()%0A%20%20%20%20%20%20%20%20n2%20%3D%20p%5B0%2C%20~positive_second%5D.sum()%0A%20%20%20%20%20%20%20%20n3%20%3D%20p%5B1%2C%20~positive_third%5D.sum()%0A%20%20%20%20%20%20%20%20scores%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20p2%20%2B%20p3%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20p2%20-%20n3%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20-n2%20%2B%20p3%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20-n2%20-%20n3%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20best_option%20%3D%20options%5Bnp.argmax(scores)%5D%0A%20%20%20%20%20%20%20%20%23%20determine%20point%20from%20cloud%20that's%20closest%20to%20the%20vertex%20we%20found%0A%20%20%20%20%20%20%20%20transformed_diag_4%20%3D%20np.array(%5B1%2C%20*best_option%5D)%0A%20%20%20%20%20%20%20%20%0A%20%20%20%20%20%20%20%20%23%20for%20topological%20reasons%20the%20new%20vertex%20must%20connect%20to%20exactly%20these%20other%20vertices%0A%20%20%20%20%20%20%20%20transformed_diagonals%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B1%2C%200%2C%200%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B0%2C%20best_option%5B0%5D%2C%200%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B0%2C%200%2C%20best_option%5B1%5D%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20transformed_diag_4%2C%0A%20%20%20%20%20%20%20%20%5D%22%22%22%0A%20%20%20%20%20%20%20%20%23%20construct%20%22cubification%22%20map%0A%20%20%20%20%20%20%20%20transformed_edge_vectors%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20d%20-%20transformed_diag_4%20for%20d%20in%20transformed_diagonals%5B%3A3%5D%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20edge_cubifier%20%3D%20np.linalg.inv(np.column_stack(transformed_edge_vectors))%0A%20%20%20%20%20%20%20%20fitted_transformation_mat%20%3D%20edge_cubifier%20%40%20estimated_transformation_mat%0A%0A%20%20%20%20%20%20%20%20cubifier%20%3D%20AffineMap(%0A%20%20%20%20%20%20%20%20%20%20%20%20fitted_transformation_mat%2C%20-edge_cubifier%20%40%20transformed_diag_4%0A%20%20%20%20%20%20%20%20).after(AffineMap(np.eye(dim_space)%2C%20-center))%0A%20%20%20%20%20%20%20%20return%20CubeFit(diagonal_orthonormalizer%2C%20cubifier)%0A%20%20%20%20return%20CubeFit%2C%20proj_cube_fit%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20An%20example%20on%20generated%20data%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20num_points%20%3D%20mo.ui.slider(50%2C%2020_000%2C%20show_value%3DTrue%2C%20value%3D10_000)%0A%20%20%20%20mo.md(f%22Number%20of%20points%20in%20generated%20sample%3A%20%7Bnum_points%7D%22)%0A%20%20%20%20return%20(num_points%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20num_points)%3A%0A%20%20%20%20n_dims%20%3D%203%0A%20%20%20%20%23%20generate%20some%20points%20in%20unit%20cube%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(1)%0A%0A%20%20%20%20%23%20solid%20cube%0A%20%20%20%20_p%20%3D%20_rng.uniform(0%2C%201%2C%20size%3D(3%2C%20num_points.value))%0A%20%20%20%20_p%20-%3D%200.5%0A%0A%20%20%20%20%23%20cube%20%22shell%22%0A%20%20%20%20%22%22%22%0A%20%20%20%20_p%20%3D%20_rng.uniform(-1%2C%201%2C%20size%3D(3%2C%2010_000))%0A%20%20%20%20_p%20%3D%20_p%20%2F%20np.linalg.norm(_p%2C%20ord%3Dnp.inf%2C%20axis%3D0)%0A%20%20%20%20_p%20%3D%200.5%20*%20_p%20%2B%200.5%0A%20%20%20%20%22%22%22%0A%0A%20%20%20%20%23%20stretch%20axes%20a%20bit%0A%20%20%20%20%23%20_p2%20%3D%20np.array(%5B2%2C%204%2C%201%5D).reshape(-1%2C%201)%20*%20_p%0A%20%20%20%20%23%20calc%20rotation%20matrix%0A%20%20%20%20_rot_ax%20%3D%20np.array(%5B1%2C%202%2C%200.7%5D)%0A%20%20%20%20_rot_ax%20%3D%20_rot_ax%20%2F%20np.linalg.norm(_rot_ax)%0A%20%20%20%20_rot_angle%20%3D%20np.pi%20*%203%20%2F%204%0A%20%20%20%20_cp_mat%20%3D%20np.array(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B0%2C%20-_rot_ax%5B2%5D%2C%20_rot_ax%5B1%5D%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B_rot_ax%5B2%5D%2C%200%2C%20-_rot_ax%5B0%5D%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5B-_rot_ax%5B1%5D%2C%20_rot_ax%5B0%5D%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20rot_mat%20%3D%20(%0A%20%20%20%20%20%20%20%20np.eye(n_dims)%0A%20%20%20%20%20%20%20%20%2B%20np.sin(_rot_angle)%20*%20_cp_mat%0A%20%20%20%20%20%20%20%20%2B%20(1%20-%20np.cos(_rot_angle))%20*%20(_cp_mat%20%40%20_cp_mat)%0A%20%20%20%20)%0A%20%20%20%20%23%20rotate%20points%20and%20skew%20axes%20a%20bit%0A%20%20%20%20transformation_mat%20%3D%20rot_mat%20%40%20np.array(%0A%20%20%20%20%20%20%20%20%5B%5B2%2C%201%2C%202.85%5D%2C%20%5B0%2C%203%2C%200.25%5D%2C%20%5B40%2C%200%2C%201.5%5D%5D%0A%20%20%20%20)%0A%20%20%20%20%23%20_center%20%3D%20np.mean(_p%2C%20axis%3D1)%0A%20%20%20%20_p%20%3D%20_p%0A%20%20%20%20points%20%3D%20(%0A%20%20%20%20%20%20%20%20transformation_mat%20%40%20_p%0A%20%20%20%20%20%20%20%20%2B%200.1%20*%20_rng.normal(size%3D_p.shape)%0A%20%20%20%20%20%20%20%20%2B%20100%20*%20_rng.normal(size%3D(3%2C%201))%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20points%20%3D%20points%20-%20np.mean(points%2C%20axis%3D1).reshape(3%2C%201)%0A%0A%20%20%20%20points_unit%20%3D%20_p%0A%20%20%20%20print(transformation_mat)%0A%20%20%20%20return%20n_dims%2C%20points%2C%20points_unit%2C%20rot_mat%2C%20transformation_mat%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Note%20that%20the%20generated%20data%20is%20quite%20heavily%20skewed%20and%20a%20bit%20noisy%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(go%2C%20points)%3A%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(dict(title%3D%22The%20generated%20data%20sample%22))%0A%0A%20%20%20%20_points%20%3D%20points%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D2%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.6%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig.update_layout(scene_aspectmode%3D%22cube%22)%0A%20%20%20%20_fig.layout.scene.camera.projection.type%20%3D%20%22orthographic%22%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(points%2C%20proj_cube_fit)%3A%0A%20%20%20%20fit%20%3D%20proj_cube_fit(points%2C%20points_centered%3DFalse)%0A%20%20%20%20return%20(fit%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(fit%2C%20mo%2C%20np%2C%20points)%3A%0A%20%20%20%20%23%20we%20want%20to%20map%20the%20points%20into%20%5B0%2C1%5D%C2%B3%20hence%20we%20should%20have%20an%20infinity-norm%20ball%20of%20radius%200.5%20around%20(1%2F2%2C%20...%2C%201%2F2).%20Points%20that%20fall%20outside%20of%20this%20ball%20are%20%22incorrect%22%0A%20%20%20%20algebraic_residual%20%3D%20np.linalg.norm(%0A%20%20%20%20%20%20%20%20fit.cubifier.apply_to(points)%20-%200.5%20*%20np.ones(3).reshape(3%2C%201)%2C%0A%20%20%20%20%20%20%20%20ord%3Dnp.inf%2C%0A%20%20%20%20%20%20%20%20axis%3D0%2C%0A%20%20%20%20)%0A%20%20%20%20res%20%3D%20np.sum(algebraic_residual%20-%200.5%2C%20where%3Dalgebraic_residual%20%3E%200.5)%0A%20%20%20%20mo.md(f%22%23%23%23%20Algebraic%20residual%20of%20fit%3A%20%7Bres%3A.3f%7D%22).center()%0A%20%20%20%20return%20algebraic_residual%2C%20res%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20**Note**%3A%20this%20residual%20doesn't%20tell%20the%20whole%20story.%20Since%20the%20distance%20function%20%24d_C%24%20of%20the%20cube%20is%20%240%24%20on%20its%20interior%20any%20large%20enough%20cube%20(that%20encompasses%20all%20of%20the%20data)%20trivially%20has%20a%20residual%20of%200.%0A%20%20%20%20%20%20%20%20So%20instead%20we%20really%20want%20to%20consider%20objectives%20of%20the%20form%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Csum_%7Bi%3D1%7D%5EN%20d_C(x_i)%5E2%20%2B%20%5Clambda%20%5Coperatorname%7Bvol%7D(C)%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20for%20regularization%20constants%20%24%5Clambda%20%3E%200%24.%20More%20on%20that%20when%20we%20get%20to%20the%20gradient%20iteration.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(fit%2C%20go%2C%20points)%3A%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(dict(title%3D%22The%20cubified%20data%22))%0A%0A%20%20%20%20_points%20%3D%20fit.cubifier%20%40%20points%0A%20%20%20%20%23%20we%20%22hollow%20out%22%20the%20cube%20a%20bit%20to%20make%20the%20plot%20more%20performant%0A%20%20%20%20%23%20_points%20%3D%20_points%5B%3A%2C%20algebraic_residual%20%3E%200.3%5D%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D2%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.6%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig.update_layout(scene_aspectmode%3D%22cube%22)%0A%20%20%20%20_fig.layout.scene.camera.projection.type%20%3D%20%22orthographic%22%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(fit%2C%20go%2C%20points)%3A%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(dict(title%3D%22Data%20with%20orthonormalized%20diagonals%22))%0A%0A%20%20%20%20_points%20%3D%20fit.diagonal_orthonormalizer%20%40%20points%0A%20%20%20%20%23%20we%20%22hollow%20out%22%20the%20cube%20a%20bit%20to%20make%20the%20plot%20more%20performant%0A%20%20%20%20%23%20_points%20%3D%20_points%5B%3A%2C%20algebraic_residual%20%3E%200.3%5D%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D2%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.6%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig.update_layout(scene_aspectmode%3D%22cube%22)%0A%20%20%20%20_fig.layout.scene.camera.projection.type%20%3D%20%22orthographic%22%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20may%20be%20able%20to%20improve%20our%20estimate%20by%20iterating%20the%20algorithm%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(fit%2C%20points%2C%20proj_cube_fit)%3A%0A%20%20%20%20estimated_map2%20%3D%20proj_cube_fit(fit.cubifier.apply_to(points)).cubifier.after(%0A%20%20%20%20%20%20%20%20fit.cubifier%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20or%20by%20utilizing%20the%20orthonormalized%20representation%0A%20%20%20%20estimated_map3%20%3D%20proj_cube_fit(%0A%20%20%20%20%20%20%20%20fit.diagonal_orthonormalizer.apply_to(points)%0A%20%20%20%20).cubifier.after(fit.diagonal_orthonormalizer)%0A%20%20%20%20return%20estimated_map2%2C%20estimated_map3%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(estimated_map2%2C%20go%2C%20points)%3A%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(dict(title%3D%22The%20doubly%20cubified%20data%22))%0A%0A%20%20%20%20_points%20%3D%20estimated_map2%20%40%20points%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D2%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.6%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig.update_layout(scene_aspectmode%3D%22cube%22)%0A%20%20%20%20_fig.layout.scene.camera.projection.type%20%3D%20%22orthographic%22%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(estimated_map2%2C%20mo%2C%20np%2C%20points)%3A%0A%20%20%20%20%23%20we%20want%20to%20map%20the%20points%20into%20%5B0%2C1%5D%C2%B3%20hence%20we%20should%20have%20an%20infinity-norm%20ball%20of%20radius%200.5%20around%20(1%2F2%2C%20...%2C%201%2F2).%20Points%20that%20fall%20outside%20of%20this%20ball%20are%20%22incorrect%22%0A%20%20%20%20_algebraic_residual%20%3D%20np.linalg.norm(%0A%20%20%20%20%20%20%20%20estimated_map2.apply_to(points)%20-%200.5%20*%20np.ones(3).reshape(3%2C%201)%2C%0A%20%20%20%20%20%20%20%20ord%3Dnp.inf%2C%0A%20%20%20%20%20%20%20%20axis%3D0%2C%0A%20%20%20%20)%0A%20%20%20%20_res%20%3D%20np.sum(_algebraic_residual%20-%200.5%2C%20where%3D_algebraic_residual%20%3E%200.5)%0A%20%20%20%20mo.md(f%22%23%23%23%20Algebraic%20residual%20of%20second%20fit%3A%20%7B_res%3A.3f%7D%22).center()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20%23%20Projected%20Gradient%20Iteration%0A%0A%20%20%20%20%20%20%20%20There's%20also%20the%20option%20of%20improving%20the%20map%20determined%20in%20the%20first%20step%20by%20a%20subsequent%20projected%20gradient%20iteration.%20However%20care%20has%20to%20be%20taken%20to%20control%20the%20step%20size.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%23%20Some%20background%0A%0A%20%20%20%20%20%20%20%20Given%20a%20convex%20set%20%24C%24%20(in%20our%20case%20%24C%20%3D%20%5B0%2C1%5D%5Ed%24%20the%20unit%20cube)%20its%20squared%20distance%20function%20%24%5Cvarphi_C(x)%20%3D%20%5Cfrac%7B1%7D%7B2%7D%20%5Cmin_%7By%20%5Cin%20C%7D%20%5ClVert%20x%20-%20y%20%5CrVert%5E2_2%24%20is%20differentiable%20and%20%24%5Cnabla%20%5Cvarphi_C(x)%20%3D%20x%20-%20P_C(x)%24%20where%20%24P_C(x)%24%20denotes%20the%20(well%20defined)%20projection%20onto%20%24C%24.%0A%0A%20%20%20%20%20%20%20%20Since%20%24C%24%20is%20a%20product%20in%20our%20case%20we%20can%20project%20componentwise%20onto%20%24%5B0%2C1%5D%24%20which%20is%20a%20simple%20clamping.%0A%0A%20%20%20%20%20%20%20%20Given%20data%20%24x_1%2C%20...%2C%20x_N%24%20we%20determine%20the%20*optimal%20cube*%20%24A%5E%7B-1%7D%20C%20%5Csubseteq%20%5CR%5E3%24%20by%20solving%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cmin_%7BA%20%5Cin%20O(3)%7D%20%5Csum_%7Bn%3D1%7D%5EN%20d_%7BA%5E%7B-1%7D%20C%7D(x_n)%5E2%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20which%20by%20the%20above%20is%20equivalent%20to%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cmin_%7BA%20%5Cin%20O(3)%7D%20%5Csum_%7Bn%3D1%7D%5EN%20%5Cvarphi_C(Ax_n)%5E2%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20which%20we%20can%20solve%20using%20projected%20gradient%20methods.%20Note%20that%20%24%5Cnabla_A%20(A%20%5Cmapsto%20%5Cvarphi_C%20(Ax_n))%7C_A%20%3D%20(Ax_n%20-%20P_C(Ax_n))%20%5Cotimes%20x_n%20%5Cin%20%5CR%5E%7B3%2C3%7D%24%20where%20%24%5Cotimes%24%20denotes%20the%20outer%20product%20of%20vectors.%0A%0A%20%20%20%20%20%20%20%20The%20orthogonal%20projection%20(w.r.t.%20Frobenius%20inner%20product)%20onto%20the%20set%20of%20orthogonal%20matrices%20is%20given%20by%20%24A%20%5Cmapsto%20UV%5ET%24%20where%20%24A%3DUSV%5ET%24%20is%20an%20SVD%20of%20%24A%24.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20This%20objective%20has%20some%20issues%20when%20we%20consider%20not%20just%20orthogonal%20matrices%3A%20as%20noted%20previously%2C%20for%20the%20general%20case%20we%20really%20want%20to%20minimize%20a%20weighted%20sum%20of%20the%20volume%20of%20the%20cube%20and%20the%20residual%20error%20of%20the%20fit.%0A%0A%20%20%20%20%20%20%20%20In%20our%20setup%20this%20translates%20to%20solving%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cmin_%7BA%7D%20%5Csum_%7Bi%3D1%7D%5EN%20d_%7B%5B0%2C1%5D%5En%7D(A%20x_i)%5E2%20%2B%20%5Clambda%20%5Coperatorname%7Bdet%7D(A).%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20i.e.%0A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cmin_%7BA%7D%20%5Csum_%7Bn%3D1%7D%5EN%20%5Cvarphi_C(Ax_n)%5E2%20%2B%20%5Clambda%20%5Coperatorname%7Bdet%7D(A).%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20As%20previously%20determined%20we%20have%20(with%20some%20abuse%20of%20notation%3B%20we're%20really%20taking%20a%20derivative%20w.r.t.%20the%20entries%20of%20the%20matrix%20here)%20%24%5Cnabla_A%20(A%20%5Cmapsto%20%5Cvarphi_C%20(Ax_n))%7C_A%20%3D%20(Ax_n%20-%20P_C(Ax_n))%20%5Cotimes%20x_n%24%20and%20by%20Jacobi's%20formula%20we%20have%20(assuming%20an%20invertible%20matrix)%20%24%5Cnabla%20%5Cdet%7C_A%20%3D%20%5Coperatorname%7Badj%7D(A)%5ET%20%3D%20%5Cdet(A)%20A%5E%7B-T%7D%24%20so%20that%20our%20full%20gradient%20becomes%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(Ax_n%20-%20P_C(Ax_n))%20%5Cotimes%20x_n%20%2B%20%5Clambda%20%5Cdet(A)%20A%5E%7B-T%7D%0A%20%20%20%20%20%20%20%20%5C%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20from%20numba%20import%20njit%0A%20%20%20%20return%20(njit%2C)%0A%0A%0A%40app.cell%0Adef%20_(njit%2C%20np)%3A%0A%20%20%20%20def%20proj_unit_cube(points%3A%20np.ndarray%2C%20out%3DNone)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20return%20np.clip(points%2C%200%2C%201%2C%20out%3DNone)%0A%0A%0A%20%20%20%20%40njit%0A%20%20%20%20def%20proj_boundary_unit_cube(points%3A%20np.ndarray%2C%20out%3DNone)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Project%20a%20set%20of%20points%20onto%20the%20boundary%20of%20the%20unit%20cube%20%5B0%2C%201%5D%5En.%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20points%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20An%20n%20x%20N%20array%20where%20n%20is%20the%20dimension%20of%20space%20and%20N%20is%20the%20number%20of%20points.%0A%20%20%20%20%20%20%20%20%20%20%20%20out%3A%20np.ndarray%2C%20optional%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20If%20provided%2C%20the%20output%20is%20stored%20in%20this%20array.%0A%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20np.ndarray%3A%20The%20projected%20points%2C%20n%20x%20N%20array.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20%23%20Clamp%20points%20to%20the%20interior%20of%20the%20unit%20cube%0A%20%20%20%20%20%20%20%20clamped%20%3D%20np.clip(points%2C%200%2C%201%2C%20out%3Dout)%0A%0A%20%20%20%20%20%20%20%20%23%20Compute%20distances%20to%20the%20boundaries%20(0%20and%201)%20for%20each%20coordinate%0A%20%20%20%20%20%20%20%20dist_to_0%20%3D%20clamped%0A%20%20%20%20%20%20%20%20dist_to_1%20%3D%201%20-%20clamped%0A%0A%20%20%20%20%20%20%20%20%23%20For%20each%20point%2C%20identify%20the%20coordinate%20closest%20to%20the%20boundary%0A%20%20%20%20%20%20%20%20min_dist_to_boundary%20%3D%20np.minimum(dist_to_0%2C%20dist_to_1)%0A%20%20%20%20%20%20%20%20idx_closest%20%3D%20np.argmin(min_dist_to_boundary%2C%20axis%3D0)%0A%0A%20%20%20%20%20%20%20%20%23%20Set%20the%20closest%20coordinate%20of%20each%20point%20to%20the%20respective%20boundary%20(0%20or%201)%0A%20%20%20%20%20%20%20%20for%20i%2C%20idx%20in%20enumerate(idx_closest)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20clamped%5Bidx%2C%20i%5D%20%3D%200%20if%20dist_to_0%5Bidx%2C%20i%5D%20%3C%3D%20dist_to_1%5Bidx%2C%20i%5D%20else%201%0A%0A%20%20%20%20%20%20%20%20return%20clamped%0A%0A%0A%20%20%20%20%40njit%0A%20%20%20%20def%20det_gradient(A)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Compute%20the%20gradient%20of%20det(A)%20with%20respect%20to%20the%20matrix%20A.%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20A%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20A%20square%20matrix%20A%20(n%20x%20n).%0A%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20np.ndarray%3A%20The%20gradient%20of%20det(A)%20with%20respect%20to%20A.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20det_A%20%3D%20np.linalg.det(A)%0A%20%20%20%20%20%20%20%20inv_A_T%20%3D%20np.linalg.inv(A).T%20%20%23%20Transpose%20of%20the%20inverse%0A%20%20%20%20%20%20%20%20return%20det_A%20*%20inv_A_T%0A%0A%0A%20%20%20%20def%20make_objective_val_and_grad(%0A%20%20%20%20%20%20%20%20points%3A%20np.ndarray%2C%20regularization_factor%3A%20float%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Points%20is%20d%20%5Ctimes%20N%20array%20of%20cartesian%20point%20coordinates%0A%20%20%20%20%20%20%20%20where%20d%20is%20the%20dimensions%20of%20space%20and%20N%20is%20the%20number%20of%20point.%0A%0A%20%20%20%20%20%20%20%20Call%20these%20x_1%2C%20...%2C%20x_N%0A%20%20%20%20%20%20%20%20%22%22%22%0A%0A%20%20%20%20%20%20%20%20def%20objective_val_and_grad(matrix%3A%20np.ndarray)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Call%20the%20matrix%20A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20This%20is%20A%20x_n%0A%20%20%20%20%20%20%20%20%20%20%20%20transformed_points%20%3D%20matrix%20%40%20points%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20proj%20%3D%20proj_boundary_unit_cube(transformed_points)%0A%20%20%20%20%20%20%20%20%20%20%20%20proj%20%3D%20proj_unit_cube(transformed_points)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20gradient%20of%20squared%20distance%20function%20at%20Ax_n%0A%20%20%20%20%20%20%20%20%20%20%20%20dist_grad%20%3D%20transformed_points%20-%20proj%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20%22gradient%22%20w.r.t%20matrix%3B%20this%20compute%20the%20outer%20products%20(Ax_n%20-%20P_C(Ax_n))%20%5Cotimes%20x_n%20%3D%20(Ax_n%20-%20P_C(Ax_n))%20x_n%5E(T)%20for%20each%20n%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20so%20this%20is%20morally%20%5Bnp.outer(dist_grad%5B%3A%2C%20i%5D%2C%20points%5B%3A%2C%20i%5D)%20for%20i%20in%20range(num_points)%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20mat_grad%20%3D%20np.einsum(%22ik%2Cjk-%3Eijk%22%2C%20dist_grad%2C%20points)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20to%20obtain%20the%20full%20gradient%20we%20simply%20sum%20over%20all%20points%0A%20%20%20%20%20%20%20%20%20%20%20%20full_grad%20%3D%20np.sum(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20mat_grad%2C%20axis%3D2%0A%20%20%20%20%20%20%20%20%20%20%20%20)%20%2B%20regularization_factor%20*%20det_gradient(matrix)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20lets%20also%20return%20the%20objective%20value%20since%20we%20can%20easily%20compute%20it%20here%0A%20%20%20%20%20%20%20%20%20%20%20%20objective_value%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20this%20should%20be%20np.linalg.vecdot(dist_grad%2C%20dist_grad%2C%20axis%3D0).sum())%20but%20can't%20use%20current%20numpy%20due%20to%20scipy%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.einsum(%22ij%2Cij-%3Ej%22%2C%20dist_grad%2C%20dist_grad).sum()%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20objective_value%2C%20full_grad%0A%0A%20%20%20%20%20%20%20%20return%20objective_val_and_grad%0A%0A%0A%20%20%20%20def%20fit_cube_naive_subgradient(%0A%20%20%20%20%20%20%20%20points%3A%20np.ndarray%2C%0A%20%20%20%20%20%20%20%20starting_guess%3A%20np.ndarray%2C%0A%20%20%20%20%20%20%20%20max_iters%3A%20int%20%3D%2010_000%2C%0A%20%20%20%20%20%20%20%20abs_tol%3A%20float%20%3D%201e-3%2C%0A%20%20%20%20%20%20%20%20step_size%3A%20float%20%3D%201e-4%2C%0A%20%20%20%20%20%20%20%20regularization_factor%3D10%2C%0A%20%20%20%20)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20%22%22%22Fits%20a%20cube%20to%20the%20given%20points%20by%20running%20a%20projected%20gradient%20method%22%22%22%0A%20%20%20%20%20%20%20%20matrix%20%3D%20starting_guess.copy()%0A%20%20%20%20%20%20%20%20obj_val_and_grad%20%3D%20make_objective_val_and_grad(%0A%20%20%20%20%20%20%20%20%20%20%20%20points%2C%20regularization_factor%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%23%20super%20stupid%20gradient%20descent%0A%20%20%20%20%20%20%20%20for%20_%20in%20range(max_iters)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20obj_val%2C%20obj_grad%20%3D%20obj_val_and_grad(matrix)%0A%20%20%20%20%20%20%20%20%20%20%20%20descent_direction%20%3D%20-obj_grad%0A%20%20%20%20%20%20%20%20%20%20%20%20grad_norm%20%3D%20np.linalg.norm(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20descent_direction%0A%20%20%20%20%20%20%20%20%20%20%20%20)%20%20%23%20frobenious%20norm%20of%20gradient%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20grad_norm%20%3C%20abs_tol%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20matrix%20%2B%3D%20step_size%20*%20descent_direction%0A%20%20%20%20%20%20%20%20%20%20%20%20svd%20%3D%20np.linalg.svd(matrix%2C%20full_matrices%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20we%20%22project%22%20(not%20really)%20onto%20the%20invertible%20matrices%20by%20forcing%20all%20singular%20values%20to%20be%20sliiiightly%20positive.%0A%20%20%20%20%20%20%20%20%20%20%20%20matrix%20%3D%20(svd.U%20*%20np.maximum(svd.S%2C%201e-4))%20%40%20svd.Vh%0A%20%20%20%20%20%20%20%20return%20matrix%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20det_gradient%2C%0A%20%20%20%20%20%20%20%20fit_cube_naive_subgradient%2C%0A%20%20%20%20%20%20%20%20make_objective_val_and_grad%2C%0A%20%20%20%20%20%20%20%20proj_boundary_unit_cube%2C%0A%20%20%20%20%20%20%20%20proj_unit_cube%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(AffineMap%2C%20estimated_map2%2C%20fit_cube_naive_subgradient%2C%20np%2C%20points)%3A%0A%20%20%20%20_mat%20%3D%20fit_cube_naive_subgradient(%0A%20%20%20%20%20%20%20%20estimated_map2.apply_to(points)%2C%0A%20%20%20%20%20%20%20%20starting_guess%3Dnp.identity(3)%2C%0A%20%20%20%20%20%20%20%20max_iters%3D50_000%2C%0A%20%20%20%20%20%20%20%20abs_tol%3D1e-4%2C%0A%20%20%20%20%20%20%20%20step_size%3D2e-3%2C%0A%20%20%20%20%20%20%20%20regularization_factor%3D10%2C%0A%20%20%20%20)%0A%20%20%20%20iterated_estimate%20%3D%20AffineMap(_mat%2C%20np.zeros(3)).after(estimated_map2)%0A%20%20%20%20return%20(iterated_estimate%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(iterated_estimate%2C%20mo%2C%20np%2C%20points)%3A%0A%20%20%20%20%23%20we%20want%20to%20map%20the%20points%20into%20%5B0%2C1%5D%C2%B3%20hence%20we%20should%20have%20an%20infinity-norm%20ball%20of%20radius%200.5%20around%20(1%2F2%2C%20...%2C%201%2F2).%20Points%20that%20fall%20outside%20of%20this%20ball%20are%20%22incorrect%22%0A%20%20%20%20_algebraic_residual%20%3D%20np.linalg.norm(%0A%20%20%20%20%20%20%20%20iterated_estimate.apply_to(points)%20-%200.5%20*%20np.ones(3).reshape(3%2C%201)%2C%0A%20%20%20%20%20%20%20%20ord%3Dnp.inf%2C%0A%20%20%20%20%20%20%20%20axis%3D0%2C%0A%20%20%20%20)%0A%20%20%20%20_res%20%3D%20np.sum(_algebraic_residual%20-%200.5%2C%20where%3D_algebraic_residual%20%3E%200.5)%0A%20%20%20%20mo.md(f%22%23%23%23%20Algebraic%20residual%20of%20projected%20gradient%20fit%3A%20%7B_res%3A.3f%7D%22).center()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(go%2C%20iterated_estimate%2C%20points)%3A%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(%0A%20%20%20%20%20%20%20%20dict(title%3D%22The%20estimate%20of%20the%20map%20after%20projected%20gradient%20iteration%22)%0A%20%20%20%20)%0A%0A%20%20%20%20_points%20%3D%20iterated_estimate%20%40%20points%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D2%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.6%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig.update_layout(scene_aspectmode%3D%22cube%22)%0A%20%20%20%20_fig.layout.scene.camera.projection.type%20%3D%20%22orthographic%22%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%20Approximating%20the%20data%0A%0A%20%20%20%20%20%20%20%20By%20inverting%20the%20determined%20maps%20we%20can%20now%20obtain%20approximations%20to%20the%20data.%0A%0A%20%20%20%20%20%20%20%20This%20meshing%20example%20assumes%20the%20%22cubifier%22%20map%20to%20be%20perfect%20which%20might%20exclude%20some%20points%20that%20are%20considered%20noise.%20To%20include%20these%20replace%20the%20cube%20vertices%20by%20the%20corresponding%20ones%20of%20the%20bounding%20box%20of%20the%20cubified%20data.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(go%2C%20iterated_estimate%2C%20np%2C%20points)%3A%0A%20%20%20%20%23%20Create%20the%203D%20scatter%20plot%0A%20%20%20%20_points%20%3D%20points%0A%20%20%20%20_fig%20%3D%20go.Figure()%0A%20%20%20%20_fig.update_layout(dict(title%3D%22Finally%20a%20mesh%20approximation%20for%20the%20data%22))%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Scatter3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_points%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_points%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_points%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mode%3D%22markers%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20marker%3Ddict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3D5%2C%20%20%23%20Size%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.8%2C%20%20%23%20Opacity%20of%20markers%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Sample%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_tets%20%3D%20np.array(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%201%2C%203)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%201%2C%205)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%202%2C%203)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%202%2C%206)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%204%2C%206)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(0%2C%204%2C%205)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(1%2C%203%2C%207)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(1%2C%205%2C%207)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(2%2C%203%2C%207)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(2%2C%206%2C%207)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(4%2C%205%2C%207)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(4%2C%206%2C%207)%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20%20%20%20%20dtype%3Dint%2C%0A%20%20%20%20).T%0A%0A%20%20%20%20_unit_cube_verts%20%3D%20%5B%5D%0A%20%20%20%20%23%20_prev_vert%20%3D%20np.array(%5B0%2C0%2C0%5D)%0A%20%20%20%20for%20_i%20in%20range(0%2C%200b111%20%2B%201)%3A%0A%20%20%20%20%20%20%20%20_unit_cube_verts.append(np.array(%5B1%20%26%20_i%2C%20(2%20%26%20_i)%20%3E%3E%201%2C%20(4%20%26%20_i)%20%3E%3E%202%5D))%0A%20%20%20%20_unit_cube_verts%20%3D%20np.array(_unit_cube_verts).T%0A%20%20%20%20_cube_verts%20%3D%20iterated_estimate.inv().apply_to(_unit_cube_verts)%0A%0A%0A%20%20%20%20_fig.add_trace(%0A%20%20%20%20%20%20%20%20go.Mesh3d(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%3D_cube_verts%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20y%3D_cube_verts%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20z%3D_cube_verts%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20i%3D_tets%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20j%3D_tets%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20k%3D_tets%5B2%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20name%3D%22Fitted%20Cube%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20opacity%3D0.5%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A